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Abstract 



In this paper, we present algorithms for computing approximate hulls and centerpoints for collections of 
matrices in positive definite space. There are many applications where the data under consideration, rather than 
being points in a Euclidean space, are positive definite (p.d.) matrices. These applications include diffusion tensor 
imaging in the brain, elasticity analysis in mechanical engineering, and the theory of kernel maps in machine 
learning. Our work centers around the notion of a horobalk the limit of a ball fixed at one point whose radius 
goes to infinity. Horoballs possess many (though not all) of the properties of halfspaces; in particular, they lack a 
strong separation theorem where two horoballs can completely partition the space. In spite of this, we show that 
we can compute an approximate "horoball hull" that strictly contains the actual convex hull. This approximate 
hull also preserves geodesic extents, which is a result of independent value: an immediate corollary is that we 
can approximately solve problems like the diameter and width in positive definite space. We also use horoballs 
to show existence of and compute approximate robust centerpoints in positive definite space, via the horoball- 
equivalent of the notion of depth. 
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1 Introduction 



There are many application areas where the basic objects of interest, rather than points in Euclidean space, are 
symmetric positive-definite n x n matrices (denoted by P(n)). In diffusion tensor imaging J3), matrices in P(3) 
model the flow of water at each voxel of a brain scan. In mechanical engineering ifTTTl . stress tensors are modeled as 
elements of P(6). Kernel matrices in machine learning are elements of P(n) ll25l . 

In all these areas, a problem of great interest is the analysis |fT3l[T4l of collections of such matrices (finding central 
points, clustering, doing regression). For all of these problems, we need the same kinds of geometric tools available 
to us in Euclidean space, including basic structures like halfspaces, convex hulls, Voronoi diagrams, various notions 
of centers, and the like. P(n) is non-Euclidean; in particular, it is negatively (and variably) curved, which poses 
fundamental problems for the design of geometric algorithms. This is in contrast to hyperbolic space (which has 
constant curvature of —1), in which many standard geometric algorithms carry over. 

In this paper, we develop a number of basic tools for manipulating positive definite space, with a focus on appli- 
cations in data analysis. 

1.1 Our Work 

Horoballs. A main technical contribution of this work is the use of horoballs as generalization of halfspaces. 
Suppose we allow a ball to grow to infinite radius while always touching a fixed point on its boundary. In Euclidean 
space, this construction yields a halfspace; in general Cartan-Hadamard manifolds (of which P(ra) is a special case), 
this construction yields a horoball. Because of the curvature of space, horoballs are not flat and the complement of 
a horoball is not a horoball. However, we show that these objects can be effectively used as proxies for halfspaces, 
allowing us to define a number of different geometric structures in P(n). 

Ball Hulls. The first structure we study is the convex hull. Apart from its importance as a fundamental primitive 
in computational geometry, the convex hull also provides a compact description of the boundary of a data set, can be 
used to define the center of a data set (via the notion of convex hull peeling depth ||23l l2l), and also captures extremal 
properties of a data set like its diameter, width and bounding volume (even in its approximate form 01). 

The convex hull of a set of points in P(n) can be naturally defined as the intersection of all convex sets containing 
the points. Alternatively, it can be defined as the set of all points that are "convex combinations" (in a geodesic 
sense) of the input points. A significant obstacle to the convex hull in P(n) is that it is not even known whether the 
convex hull of a finite collection of points in P(ra) can be represented finitely [4J. 

Another approach to defining the convex hull is via halfspaces: we can define the convex hull in Euclidean space 
as the intersection of all halfspaces that contain all the points. Unfortunately, even this notion fails to generalize: the 
relevant structures are called totally geodesic submanifolds, and we cannot guarantee that any set of d + 1 points 
admits such a submanifold passing through them. 

Our main technical contribution here is a generalization of the convex hull called the ball hull that is based on 
the relationship between horoballs and halfplanes. The ball hull is the intersection of all horoballs that contain the 
input points. Although the ball hull itself might require an infinite number of balls to define it, it is closed, it can be 
approximated efficiently, it is identical to the convex hull in Euclidean space, and it always contains the convex hull 
in P(n). In the process of proving this result, we also develop a generalized notion of extent [lj in positive definite 
space that might be of independent interest for other analysis problems. 

Centerpoints. One important motivation for studying collections of points in positive definite space is to compute 
measures of centrality (or mean shapes) 1 14]. A robust centerpoint can be obtained by finding a point of maximum 
(halfspace) depth among a collection of points. We first prove, using a generalization of Helly 's theorem to negatively 
curved spaces, that for any set of points in P(n) , there exists a point of large depth, where depth is defined in terms of 
horoballs. We then develop an algorithm to compute an approximation to such a point, using an LP-type framework. 
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The point we compute is a geometric approximation: it does not approximate the depth of the optimal point, but is 
guaranteed to be close to such a point. 

1.2 Related Work 

The mathematics of Riemannian manifolds, Cartan-Hadamard manifolds and P(n) is well-understood: the book 
by Bridson and Haefliger (6l is an invaluable reference on metric spaces of nonpositive curvature, and Bhatia (U 
provides a detailed study of P(n) in particular. However, there are many fewer algorithmic results for problems in 
these spaces. To the best of our knowledge, the only prior work on algorithms for positive definite space are the 
work by Moakher II2T1I on mean shapes in positive definite space, and papers by Fletcher and Joshi fT3Tl on doing 
principal geodesic analysis in symmetric spaces, and the robust median algorithms of Fletcher et al lfl4ll for general 
manifolds (including P(n) and SO(n)). 

Geometric algorithms in hyperbolic space are much more tractable. The Poincare and Klein models of hyperbolic 
space preserve different properties of Euclidean space, and many algorithm carry over directly with no modifica- 
tions. Leibon and Letscher [18] were the first to study basic geometric primitives in general Riemannian manifolds, 
constructing Voronoi diagrams and Delaunay triangulations for sufficiently dense point sets in these spaces. Epp- 
stein [12J described hierarchical clustering algorithms in hyperbolic space. Krauthgamer and Lee lfT6l studied the 
nearest neighbor problem for points in <5-hyperbolic space; these spaces are a combinatorial generalization of neg- 
atively curved space and are characterized by global, rather than local, definitions of curvature. Chepoi ef (HI |9l 
advanced this line of research, providing algorithms for computing the diameter and minimum enclosing ball of 
collections of points in 5-hyperbolic space. 

2 Preliminaries 

P(n) is the set of symmetric positive-definite real matrices. It is a Riemannian metric space with tangent space at 
point p equal to S(n), the vector space of symmetric matrices with inner product {A, B) = ti(p~ 1 Ap~ l B). The 
exp map, exp : S(n) — * P(n) is defined exp (L4) = c(t) = pe tpA , where c(t) is the geodesic with unit tangent A 
and c(0) = p. For simplicity, we often assume that p = I so expj(tA) = e tA . The log map, log„ : P(n) — ► S(n), 
indicates direction and distance and is the inverse of exp p . The metric d(p, q) = \\ logJq) \\ = y / tr(log(p~ 1 q) 2 ). 

Convex Hulls in P(n). P(n) is an example of a proper CAT(O) space [6, 11.10], and as such admits a well-defined 
notion of convexity, in which metric balls are convex. We can define the convex hull C(X) of a set of points X as 
the smallest convex set that contains the points. This hull can be realized as the limit of an iterative procedure where 
we draw all geodesies between data points, add all the new points to the set, and repeat. 

Lemma 2.1 (0). IfX = X and X l+1 = \J afieXi [a, b], then C{X) = \JZ X { . 

Proof. We will use the notation Xoa = |J°^ Xi. It is easy to demonstrate by straightforward induction that X^ is 
contained in any convex set that contains X. Therefore C{X) D X^. 

We also know that if p, q G X^ there must be some m for which p, q G X m , since Xoo is the nested union of Xj. 
Then [p, q] C X m+ i C X^. This means that X^ is convex, so C(X) C X^. □ 

Berger [4 ] notes that it is unknown whether the convex hull of three points is in general closed, and the standing 
conjecture is that it is not. The above lemma bears this out, as it is an infinite union of closed sets, which in general 
is not closed. These facts present a significant barrier to the computation of convex hulls on general manifolds. 

2.1 Busemann Functions 

In R d , the convex hull of a finite set can be described by a finite number of hyperplanes each supported by d points 
from the set. A hyperplane through a point may also be thought of as the limiting case of a sphere whose center has 
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been moved away to infinity while a point at its surface remains fixed. We generalize this notion with the definition 
of a Busemann function. 

For this notion to work, we must restrict ourselves to a class of spaces called CAT(O) spaces. They are metric 
spaces with non-positive curvature. Additionally, they must be complete; that is, Cauchy sequences in the space 
must converge to a point in the space. Euclidean space, hyperbolic space, and P(n) are all examples of complete 
CAT(O) spaces. To talk about "sending a point away to infinity," we must provide a rigorous definition of what we 
mean by infinity in a complete CAT(O) space. 

Two geodesic rays c\, C2 : M + —>■ M in a complete CAT(O) space M are asymptotic if lim^oo d(ci(t), C2(t)) < 
K for some K G M + . If c\ and C2 are asymptotic, then we say c\ ~ C2. This forms an equivalence relation ~ so 
that [c] describes the set of all geodesies d such that c ~ c' . Let £ = [c] where £ is identified with the limit of any 
geodesic ray asymptotic to c. We say that £ is a point at infinity. Moreover, for any point x G M we can find a 
member of [c] that issues from x [6, II.8]. 

Definition 2.1. For a complete CAT(0) space M, given a geodesic ray c(t) : M + — ► M, a Busemann function 

b c : M — > R is defined 

b c (p) = lim d(p,c(t)) - t. 

t— >oo 

It should be noted that if we construct a Busemann function from any geodesic ray in [c] , it is the same function up 
to addition by a constant ||6j II.8]. It's convenient then to normalize a Busemann function by assuring that b c (I) = 0. 

A Busemann function is an example of a horofunction (H II.8]. A horosphere S r (h) C M is a level set of a 
horofunction h; that is, S r (h) = /i _1 (r), where r G M. A horoball B r (h) C M is a sublevel set of h; that is, 
B r (h) = oo, r]). Horofunctions are convex (6l II.8], so any sublevel set of a horofunction is convex, and 

therefore any horoball is convex. 

Example: Busemann functions in M n . As an illustration, we can easily compute the Busemann function in 
Euclidean space associated with a ray c(t) = tu, where u is a unit vector. Since Hindoo <fe(||p — £u|| + t) = 1, 

b c (p) = lim (\\p - tu\\ - t) 

i— >oo 

= lim ±-(\\p-tu\\ 2 -t 2 ) 

t^oo It 

= lim ^(|H| 2 -2(p,iu) + ||tu|| 2 -i 2 ) 
\\p\\ 2 

= } im -^T ~ VP> u ) = - (P. u ) • 

t^oo 2r 

Horospheres in Euclidean space are then just hyperplanes, and horoballs are halfspaces. 
2.1.1 Decomposing P(n) 

In order to construct Busemann functions in P(n) it is necessary to decompose the space into simpler components. 
The notion of a horospherical projection will be very useful. 

The horospherical group. There is a subgroup of GL(ra), (the horospherical group), that leaves the Buse- 
mann function b c invariant [6, 11.10]. That is, given p G P(n), and v G Ng, b c {vpv T ) = b c (p). Let A be diagonal, 
where An > Ajj, Mi > j. Let c(t) = e tA , and £ = c(oo). Then ^ G if and only if ^ is a upper- triangular matrix 
with ones on the diagona Q If A G S(n) is not sorted-diagonal, we may still use this characterization of N% without 
loss of generality, since we may compute an appropriate diagonalization A = QA'Q T , QQ T = I, then apply the 
isometry Q T pQ to any element p G P(n). 

'For simplicity, we consider only those rays with unique diagonal entries, but this definition may be extended to those with multiplicity. 
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Figure 1: Left: projection of X C P(2) onto det(ir) = 1. Right: X C P(2). Two horospheres are drawn in both views. 

Flats. Let A G S(n) and c(t) = e tA as above. If we consider all elements / G P(n) that share eigenvectors Q with 
e A , then all such elements commute with each other and fe A = e A f. We call this space F, the n-flat containing c. 
Since we may assume that Q € SO(n), every flat F corresponds to an element of SO(n). Moreover, since members 
of F commute, log(ah) = log a + log b for all a,b G F. So if u and v are in F, then the distance between them is 
y / tr(log(u~ 1 v) 2 ) = -y/tr((log(t;) — log(u)) 2 ). Since y / tr((-) 2 ) is a Euclidean norm on log(F), we have that F is 
isometric to W 1 with a Euclidean metric under log(-). 



Horospherical projection. Given p G P(n), there is a unique decomposition p = vfv T where (v, f) G N^xF |6j 
11.10]. Let p G P(n) and (z/, /) G x F. If p = vfv T , then define the horospherical projection function 
TTp : P(n) — > .F as ttf{p) = v~ l pv~ T = /. 



2.1.2 Busemann functions in P(ra). 

We can now give an explicit expression for a Busemann function in P(n). For geodesic c(t) = e , where A G S(n), 
the Busemann function b c : P(n) — > M is 

6 c (p) = -tr(Alog(7r F (p))), 

where 7ri7 is defined as above |[6j 11.10]. 

In P(2) it is convenient to visualize Busemann functions through horospheres. We can embed P(2) in R 3 where 
the log of the determinant of elements grows along one axis. The orthogonal planes contain a model of hyperbolic 
space called the Poincare disk that is modeled as a unit disk, with boundary at infinity represented by the unit circle. 
Thus the entire space can be seen as a cylinder, as shown in Figure [T] Within each cross section with constant 
determinant, the horoballs are disks tangent to the boundary at infinity. 



3 Ball Hulls 



We now introduce our variant of the convex hull in P(n), which we call the ball hull. For a subset X C P(n), the 
ball hull B(X) is the intersection of all horoballs that also contain X: 

B{X) = p| B r (b c ), XcB r (b c ). 

b c ,r 

Note that the ball hull can be seen as an alternate generalization of the Euclidean convex hull (i.e. via intersection of 
halfspaces) to P(n). Furthermore, since it is the intersection of closed sets, it is itself guaranteed to be closed. 



3.1 Properties Of The Ball Hull 

We know that any horoball is convex. Because the ball hull is the intersection of convex sets, it is itself convex 
(and therefore C(X) C B(X)). We can also show that it shares critical parts of its boundary with the convex hull 
(Theorem |3.1|>, but unfortunately, we cannot represent it as a finite intersection of horoballs (Theorem|3.2[). 
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Theorem 3.1. Every x G X (X finite) on the boundary ofB(X) is also on the boundary ofC(X) (i.e., XDdB(X) C 

xndC{x)). 

Proof. Since X C C(X), either x £ dC(X) or x £ int C(X). Assume that x £ int C(X). Then there is a 
neighborhood U of x contained wholly in C(X). Because x E c?<B(X), there is a horof unction /i such that X C 
B r (h) and = r. Since B r (h) is convex, [/ C C(X) C B r (h). This implies that h(x) < r, but = r, a 
contradiction. Thus x G 9C(X). □ 

Theorem 3.2. /« general, the ball hull cannot be described as the intersection of a finite set of horoballs. 

Proof. We construct an example in P(4) with a point set X = {xi,X2} of size 2 where the ball hull cannot be 
described as the intersection of a finite number of horoballs. Let X C H 3 , the three dimensional hyperbolic space, 
as embedded in P(4). In particular, let the geodesic that contains x\ and X2 also contain /, the identity matrix, at 
their midpoint on the geodesic. 

Consider the family of horofunctions such that for h £ "K, h(x±) = hixz). By construction, h(x\) = h{x2) = 
r for some r £ M. In the Poincare ball model of H 3 , the horospheres S r (h) are literally spheres that are tangent 
to the boundary at infinity, and touch x\ and X2- So constructed, any horosphere in !H will contact the boundary at 
infinity on a great circle that is equidistant from both points. 

The ball hull B{X) is defined {(~) B r (h) \ h £ %}, and is a "spindle" (a three dimensional lune) with tips at x\ 
and X2 and bulges out from the geodesic segment between them. Any finite family of horofunctions will intersect 
in a region strictly larger than B(X), so every horoball generated by a member in IK is necessary for B(X), and so 
there is no finite set of horoballs that describe B(X). □ 



4 The £-Ball Hul 



Theorem 3.2 indicates that we cannot maintain a finite representation of a ball hull. However, as we shall show in 



this section, we can maintain a finite-sized approximation to the ball hull. Our approximation will be in terms of 
extents: intuitively, we say that a set of horoballs approximates the ball hull if a geodesic traveling in any direction 
traverses approximately the same distance inside the ball hull as it does inside the approximate hull. 

In Euclidean space, we can capture extent by measuring the distance between two parallel hy- 
perplanes that sandwich the set. Since measuring extent by recording the distance between two 
parallel planes does not have a direct analogue in P(n), we define a notion we call a horoextent. Let 
c(t) = qe tq be a geodesic, and X C P(n). The horoextent E C (X) with respect to c is defined as: 



E C {X) 



max b c+ (p) + max 6 C _ (p) 




horoextent 



where b c + is the Busemann function created when we allow t to approach positive infinity as normal, while fe c _ is 
the Busemann function created when we allow the limit to go the other direction; that is: 



\p) = lim (d(c(t),p) 

t— >+oc 



*)> &c-(p) 



lim (d(c(t),p) +t). 



Observe that for any c, E C (X) = E C (C(X)) = E C (B(X)). (Note that we cannot simply substitute min6 c+ for 
max 6 C _ ; the horoballs are generated by Busemann functions tied to opposite points at infinity.) 

If we were to compute E C (X) in a Euclidean space, it would be apparent that the extent would be the distance 
between two parallel planes. Because the minimum distance between two Euclidean horoballs is a constant, no 
matter where we place the base point q, we can measure horoextent simply by measuring width in a particular 
direction. However, this is not true in general. For instance in P(n), horofunctions are nonlinear, so the distance 
between opposing horoballs is not constant. The width of the intersection of the opposing horoballs is taken along 
the geodesic c, and a geodesic is described by a point q and a direction A. We fix the point q so that we need only 
choose a uniform grid of directions A for our approximation. 
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Definition 4.1. An intersection of horoballs is called an e-ball hull with origin q (B £A (X) ) if for all geodesic rays c 
such that c(0) = q, \E c (B e ^ q (X)) — E C (X)\ < e. When q is clear, we will refer to B £i q(X) as just an e-ball hull. 

Shifting the origin. Let the geodesic anisotropy lf2TTl of a point p G P(n) be defined as GA(p) = d( y / det(p)I, p) 
(so if det(p) = 1 then GA(p) = d(I,p)). Let dx = max pg x d(p, I) > max pg j GA(p). The size of the e-ball 
hulls we construct will depend dx, but this is not an intrinsic parameter of the data, since we can change it merely 
by isometricaily translating the point set. For some point q G C(X), if we could translate the data set so that q was 
at the origin /, then dx < diam(X) = maxp )?g x d(p, q). We now prove that such a translation is always possible. 

Lemma 4.1. For a point q G P{n), a geodesic c such that c(0) = q, and a point set X, if we define an operation S 
on a set S £ P(n) such that p = q~ 2pq~ 5 for any p G S, then 

E C (X)=E £ (X). 

Proof Let c(t) = qe tq ' lA . Then c(t) = q~^ (qe tq ~ lA )q-^ = g5e**~ lj4 g~§ = e *«~ 2A «~ 2 = e tA . Note that 
\\A\\q = yj tr ((q~ 2 Aq~ 2 ) 2 ) = \\A\\j, so c is a geodesic such that c(0) = I with the same speed as c. Let b c be the 
Busemann function of c, so b c (p) = lim t ^ 00 (d(c(t),p) — t). Since conjugation by q~2 is an isometry on P(n), 

b c (p) = lim (d(c(t),p) -t)= lim (d{c{t),p) -t) = b t {p), 

t—>oo t—>oo 

and therefore 

E C (X) = E 8 (X). 

□ 

For convenience, we will now assume that our data has been shifted into a reasonable frame where some point 
q G C(X) is the base point of our horofunction. That is in this shifted frame / G C(X) and, we can bound, 
d x < diam(X). 



Main result. The main result of this section is a construction of a finite-sized e-ball hull. 

Theorem 4.1. ForasetX C P{n) of size N (for constant n), we can construct an e-ball hull of size 0((sinh(dx) / 
M™/ 2 J) in time 0((sinh(d x )/e)™~ 1 ( ArLn/2j + NlogN)). 



Proof Overview. We make much use of the structure of flats in our proof, so it is helpful to describe some 
conventions. Consider the set of unit-length tangent vectors at /, part of the tangent space S(n); in other words, the 
set of "directions" from /. If we choose a flat F to work in, then the tangent space of F contains a subset of those 
directions. All these directions, though, share the rotation Q identified with F. So in much of the rest of the paper, 
we refer to this rotation Q as a "direction," even though it is not a member of the tangent space. 

Our proof uses two key ideas. First, we show that within a flat F (i.e., given a direction Q G SO(n)) we can 
find a finite set of minimal horoballs exactly. This is done by showing an equivalence between halfspaces in F and 
horoballs in P(n) in Section 4. 1 The result implies that computing minimal horoballs with respect to a direction Q 
is equivalent to computing a convex hull in Euclidean space. 

Second, we show that instead of searching over the entire space of directions SO(n), we can discretize it into 
a finite set of directions such that when we calculate the horoballs with respect to each of these directions, the 
horoextents of the resulting e-ball hull are not too far from those of the ball hull. In order to do this, we prove a 
Lipschitz bound for horofunctions (and hence horoextents) on the space of directions. Since any two flats F and F' 
are identified with rotations Q and Q' , we can move a point from F to F' simply by applying the rotation Q T Q', 
and measure the angle 9 between the flats. If we consider a geodesic c C F such that c(0) = /, we can apply Q T Q' 
to c to get c', then for any point p G P(n) we bound | b c (p) — b c i (p) \ as a function of 9. 
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Proving this theorem is quite technical. We first prove a Lipschitz bound in P(2), where the space of directions is 
a circle (as in the left part of Figure [TJ. After providing a bound in P(2) we decompose the distance between two 
directions in SO(n) into (") angles defined by 2 x 2 submatrices in an n x n matrix. In this setting it is possible 
to apply the P(2) Lipschitz bound times to get the full bound. The proof for P(2) is presented in Section 



and the generalization to P(n) is presented in Section 4.3 Finally, we combine these results in an algorithm in 
Section 03] 



4.2 



The following lemma describes how geodesies (and horofunctions) are transformed by a rotation. 
Lemma 4.2. For a point p 6 P(n), a rotation matrix Q, geodesies c(t) = e tA and c'(t) = e l Q A Q T t 

b c >(p) = b c (Q T P Q). 
Proof. If A' = QAQ T is the tangent vector of c', and F' is the flat containing cf, 

b c ,{p) = -tr(^log(7r^(p))) = -ti{{QAQ T )\og{{Qv- l Q T )p{Qv" l Q T ) T )) 

= -ti{AQ T \og{Qu~\Q T pQ)u~ T Q T )Q) = -ti{A\og{u-\Q T pQ)u~ T )) 
= b c {Q T pQ). 

□ 

In particular, this allows us to pick a flat where computation of b c is convenient, and rotate the point set by Q to 
compute b c instead of attempting computation of b c > directly, which may be more cumbersome; we will utilize this 
idea later. 



4.1 Projection to A>flat 

For the first part of our proof, we establish an equivalence between horospheres and halfspaces. That is, after we 
compute the projection of our point set, we can say that the point set lies inside a horoball B r (b c ) if and only if its 
projection lies inside a halfspace H r of F (recall that F is isometric to a Euclidean space under log). 

Lemma 4.3. For any horoball B r (b c ), there is a halfspace H r C log(i ? ) C S(n) such that log(TTp(B r (b c ))) = H r . 

Proof. If b c (p) < r, p G P(n), and c(t) = e tA , then — ti(Alog(ir f{p))) < r. Since ttf{p) is positive-definite, 
log(7Tp(p)) is symmetric. But tr((-)(-)) defines an inner product on the Euclidean space of symmetric n x n 
matrices. Then the set of all Y such that — tr(AY) < r defines a halfspace whose boundary is perpendicular to 
A. ' □ 

This gives us a means to compute horoballs by using ttf to project our point set onto F, and leverage a familiar 
Euclidean environment. 

4.2 A Lipschitz bound in P(2) 
4.2.1 Rotations in P(2) 

We start with some technical lemmas that describe the locus of rotating points in P(2). 

Lemma 4.4. Given a rotation matrix Q £ SO{2) corresponding to an angle of 9/2, Q acts on a point p £ -P(2) via 
QpQ T as a rotation by 6 about the (geodesic) axis e t! = e t L 

Proof. If p = e l l, J 6 I, then QpQ T = e t QIQ T = e f I, so e l I is invariant under the action of Q. Any action 
GpG T where G £ GL(ra) is an isometry on P(n), so the distance from p to the axis e*I remains fixed [6, 11.10]. 
Computing QpQ T as a function of 9, we get: 

V ^ sin0 + ™ cos0 ^ - ^cos6 + wsm9J ' 
which is 2tt -periodic. (In fact, it is easy to see a rotation in the "coordinates" and w.) □ 
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By Lemma |4.4| we know that as we apply a rotation to p, it moves in a circle. Because any rotation Q has 
determinant 1, det(QpQ T ) = det(p). This leads to the following corollary: 

Corollary 4.1. In P{2), the radius of the circle that p travels on is GA(p) . Such a circle lies entirely within a 
submanifold of constant determinant. 

In fact, any submanifold P(2) r of points with determinant equal to some r E M + is isometric to any other such 
submanifold P(2) s for s G M + . This is seen very easily by considering the distance function ti(log(p~ 1 q)) — the 
determinants of p and q will cancel. 

We pick a natural representative of these submanifolds, P(2)i. This submanifold forms a complete metric space 
of its own that has special structure: 



Lemma 4.5. P(2)i has constant sectional curvature 



2- 



J , det(p) = xy - w 2 = 1. Let u = S±2 and v = ^ so that x = u + v 

and y = u — v. Then det(p) = u 2 — v 2 — w 2 = 1. This describes a hyperboloid of two sheets, and restricting u > 0, 
is a model for hyperbolic space H 2 . 

To analyze the metrics between the two spaces, we may consider our other point to be the identity matrix, since 
in P(2)i, d(p, q) = d(q~ 1 ^ 2 pq~ 1 ^ 2 , 1). The equivalent point on the hyperboloid to I is (n, v, w) = (1, 0, 0). The 
distance between the two points in the hyperbolic metric is 



djp(p, I) = cosh 1 (uiU2 — V\V2 — W1W2) = cosh 1 f 
And in the metric of P(n) : 



d m (p, I) = ^tr(log 2 (p)) = ^ln 2 A 1+ ln 2 A 2 = ^2 In X 1 = V2 In + ^(*±^j - 

Since this is a constant multiple of the hyperbolic metric, P(2)i is a complete metric space of constant negative 
sectional curvature. We can find the curvature k by solving 1 / ^J— k = y/2 to get k = —1/2 [6 , 1.2]. □ 



4.2.2 Bounding ||V6 C || 

To bound the error incurred by discretizing the space of directions, we need to understand the behavior of b c as a 
function of a rotation Q. We show that the derivative of a geodesic is constant on P(n). 

Lemma 4.6. For a geodesic ray c(t) = e tA , \\A\\ = 1, then ||V6 C || = 1 at any point p G P(n). 

Proof. Let u = b c (p). If we define the map 7 P to map p to its projection onto any horoball £> M _ t (6 c ), where t > 0, 
then 7 P is a geodesic ray @ II.8]. Since any j p (t) is the projection onto the horoball B u _ t (b c ), the geodesic segment 
[p> 7p(*)] i s perpendicular to B u ^ t {b c ) at 7 p (t), and therefore the tangent vector of 7 P points directly opposite to Vb c 
at 7 p (i). Because 7 P intersects B u _ t (bc) at t, b c {^ p {t)) must change at a rate opposite to 7 P (i), along j p , and since 
V6 C points in the opposite direction as 7 p (i) at 7 p (t), V6 C = — 7«(t). 

Also, since d(p, B u _ t {b c )) = t for any p G B u (b c ), for any n, and ||7 p || is constant along j p , \\Vb c \\ is the same 
anywhere in P(n). Since c(t) is the projection of c(0) onto B-t{b c ) by construction, and geodesies are unique in 
P(n), ||V6 C || = \\A\\ = 1. 

□ 
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4.2.3 A Lipschitz condition on Busemann functions in P(2) 



We are now ready to prove the main Lipschitz result in P(2). We start with a more specific bound that depends on 
the geodesic anisotropy of a point: 



Lemma 4.7. Given a point p G P{2), a rotation matrix Q corresponding to an angle of 9/2, geodesies c(i) — < 
and c'(t) = e «^Q T , 

{GA(p)\ 



tA 



Mp)-b d {p)\ < |^| - V2sinh 



V V2 J' 



Proof. The derivative of a function / along a curve ^(t) has the form (V/| 7 m, 7'(i)), and has greatest magnitude 
when the tangent vector ^'(t) to the curve and the gradient V/| 7 ( t ) are parallel. When this happens, the derivative 
reaches its maximum at ||V/| 7 (t)|| • ||7 ; (i)||. 

Since ||V6 C || = 1 anywhere by Lemma 4.6 the derivative of b c along 7 at j(t) is bounded by ||Y(£)||. We are 
interested in the case where 7(0) is the circle in P(2) defined by tracing Q(9/2)pQ(9/2) T for all — ir < 9 < n. 
By Corollary 4.1 we know that this circle has radius GA(p) and lies entirely within a submanifold of constant 
determinant, which by Lemma 4.5 also has constant curvature k = — 1/2. 

This implies that 

||V(0)|| = -^L=sinh(^r) = \/2sinh f^=p) ■ 
for any value of 9 G (— tt, n] \6, 1.6]. Then 

\bc(p) ~ b c '(p)\ = MP) - b c {Q T pQ)\ < \9\ ■ V2 S inh f^^J • 

□ 

We can now state our main Lipschitz result in P(2). 
Theorem 4.2. For Q G 50(2) corresponding to 9/2, and let j(t) = e tA and ^'(t) = e tQTAQ . Then for any p G X 

d x 



|6 7 (p) - by(p)\ < \9\ ■ \/2smh 



V2 



4.3 Generalizing to P(n) 

Now to generalize to P(n) we need to decompose the projection operation ttf{-) and the rotation matrix Q. We can 
compute 7tp recursively, and it turns out that this fact helps us to break down the analysis of rotations. Since we can 
decompose any rotation into a series of 2 x 2 rotation matrices, decomposing the computation of ttf in a similar 
manner lets us build a Lipschitz condition for P(n). 

Lemma 4.8. Let F C P (n) be the flat containing diagonal matrices, let r + s = n, and let ixp^ r and ttf,s be the 
projection operations for r x r and s x s flats of diagonal matrices, respectively. Then 

V Kf,s{Ps)J 
where p s is the lower right s x s block ofp, and h r is the Schur complement ofp s . 

Proof. TTpis computed by decomposing p into v~ 1 fi>~ T , where / is diagonal (and positive-definite) and v is upper- 
triangular, with ones on the diagonal. This can be done by computing the Schur complement of some lower right 
square block of p, and putting the complement in the upper right block of a new matrix, with the lower block in the 
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other corner. The original matrix can be reconstructed by conjugating by an upper- triangular matrix of appropriate 
form: 

' p r a\ (l r ap" 1 ^ ( h r 



P ~ \a T PsJ~\ Is A V.) Wa T h 

Performing this process recursively yields a product of upper-triangular matrices with ones on the diagonal, which 
is again an upper-triangular matrix with ones on the diagonal, yielding v, and a diagonal matrix, / = np(p). □ 

If we wish to compute ttf{p) for other flats, then we can apply the rotation of F to p, compute using the recursive 
formula described above, and then apply the opposite rotation to the resulting diagonal matrix. Most of the time, 
however, it is most convenient to condition our data so that 7Tp is computed for a diagonal flat F. 

We now wish to analyze a simpler form of rotation, one that can be broken into rotations on separate axes. 



Lemma 4.9. Given a point p 6 P{n), a rotation matrix Q = ( ^ r Q ) , suc ^ t ^ mt r + s = n and Q r , Q s are 



r x r 



Q s 

; s x s rotation matrices, respectively, and a geodesic c(t) = ef-Q A Q T with A = ^ \ sorted-diagonal, 



Hp) = -tr(A r log(7T^ r (Q^/i r Q r ))) -tx(A s log(7r F , s (Q^PsQs))) 



Proof. From Lemma 4.2 we know that b c (p) = — tx(A log(7rp((3 pQ)))- Compute ttf{Q pQ) by first decompos- 
er a \ 
a T PsJ' 

Q T pQ 



Qj \ ( Pr a\ (Q r \ _ ( QjprQr QjdQs 



p r a 
mgp^ I T 



Qs J \ a PsJ V Qs) \Qs a Qr Qs PsQ 

Now compute the Schur complement of Qjp a Q B : 

QjprQr - Qj.aQ s {Q T s p s Q s )~ l Q T s a T Q r = Q T r PrQr - Q^ aQ s Q T g p~ l Q S Q T S a T Q r 

= Q r PrQr ~ Q r a Q 

= Qr(Pr ~ ap~ 1 a T )Q r . 
But h r = p r — ap~ l a T is just the Schur complement of p s , so 

bM = -tr(Alo g (MQ T pQ))) = -tr( Arl ° g ^ A ^ fri r 



A s log(7TF, s (Qip s Q s )) y 
tv(A r log{Tr F;r {Qjh r Q r ))) - tr(A s log(ir F>s (QTp s Q s ))). 



□ 



This allows us to break the Lipschitz bound into smaller pieces that we can analyze individually. The following 
two corollaries give us a way to analyze the effects of 2 x 2 rotation matrices: 

//, \ 

Corollary 4.2. Given a point p G P(n), a rotation matrix Q = I Q I , where r + s + 2 = n, Q' is a 2 x 2 

\ ' Is) 

rotation matrix corresponding to an angle of 9/2, geodesies c(t) = e tA and c'(t) = e t( ^ A< ^ , then \b c (p) — b c '(p)\ is 
bounded as in Lemma\4.7\ 



Proof. This is easily seen after observing that I r and I s are also rotation matrices, so Lemma 4.9 can be applied 
twice. □ 
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Corollary 4.3. The results of Corollary 4.2 extend to rotations between any two coordinates, that is, where Q' is of 
the form 

'cos(6>/2) -sin(6>/2)' 
h 

v sin(fl/2) cos(#/2) 

Proof. First observe that a rotation matrix Qij that rotates between axes i and j is equal to a matrix Ei + i^Qi^ + \Ef +1 ■, 
where Eij is a permutation that moves row i to row j and shifts the intervening rows up. We assume E = -Ei+ij, 
Q' = Qi,i+\, an d Q = EQ'E T from here on. 

Assuming that A is sorted-diagonal, we can compute 6 c '(p) as: 

bAp) = -tr((EQ'E T )A(EQ'E T ) T log(((EQ'E T )v-HEQ'E T ) T )p^ 
= - tr((E T AE) log((E T u- 1 E)Q' T (E T P E)Q'(E T u~ 1 E) T )) 
= -ti(Alog(i>- 1 Q' T pQ'i>- T )) 
= -tr(ilog(7r^(Q /T pQ'))), 

which can be computed as above; some care must be taken, however, since the order of elements of A is different 
than that of A. That is, in certain places, the Schur complement of the upper corner must be taken to compute irp, 
rather than that of the lower corner. □ 

4.3.1 A Lipschitz condition on Busemann functions in P(n) 

We are now ready to prove the main Lipschitz result in P(n). We start with a more specific bound that depends on 
the distance from a point p to /: 



Lemma 4.10. Given a point p £ P{n), a rotation matrix Q E SO{n) corresponding to an angle of 9/2, geodesies 
c{t) = e tA and c'(t) = e t( ^ A ^ ', 



Mp)-bAp)\<\o\- 



Proof. Every rotation Q may be decomposed into a product of rotations Q = Q1Q2 ■ ■ - Qk where k = (™) and Qi 
is a 2 x 2 sub-block rotation corresponding to an angle of 6i/2 with \9i\ < \9\. Then 



\b c (p) ~ b c '(p)\ 



i=i 



<Ei£" l (p)-£(p)i. 



i=l 



where b a c ,(p) = b c (p) and b l c ,(p) is b c (p) with the first i rotations successively applied, so if Q[ = YYj=i Qj> 

b i A P ) = b c {{Q' l ) T p{Q' l )). 

But then 



\b i c 7 1 {p)-bl / {p)\ < \0i\ -V2smh( 



dip, I) 



and therefore 



MP) - M*)l < (t • V2sinh (^D-) < |»| ■ Q • V2 S inh 
since for all % we have \9;\ < \9\. □ 
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We can now state our main Lipschitz result in P(n). 



Theorem 4.3 (Lipschitz condition on Busemann functions in P(n)). Consider a set X C P{n), a rotation matrix 
Q G SO(n) corresponding to an angle 6/2, geodesies c(t) = e tA and c'{t) = e tc ^ A( ^ . Then for any p G X 



\b c (p)-bA P )\<\0\-(™yV2sinh(^\ 



4.4 Algorithm 



For X C P(n) we can construct e-ball hull as follows. We place a grid G £ on SO(n) so that for any Q' G SO(n), 
there is another Q G G £ such that the angle between Q and Q' is at most (e/2)/(2(2)\/2sinh((ix/v / 2))- For each 
Q G G £ , we consider ttf(X), the projection of X into the associated n-flat F associated with Q. Within F, we 
construct a convex hull of ttf(X), and return the horoball associated with each hyperplane passing through each 



facet of the convex hull, as in Lemma 4.3 



To analyze this algorithm we can now consider any direction Q' G SO(n) and a horofunction b c > that lies in 
the associated flat F'. There must be another direction Q G G £ such that the angle between Q and Q' is at most 
(e/2) j (2 (™) y/2 sinh(<ix / \/2))- Let b c be the similar horofunction to b c >, except it lies in the flat F associated with 
Q. This ensures that for any point p G X, we have \b c i{p) — b c (p)\ < e/2. Since E c i(X) depends on two points in 
X, and each point changes at most e/2 from b c > to b c we can argue that \E C /(X) — E C (X)\ < e. Since this holds 
for any direction Q' G SO(n) and for Q G G £ the function E C (X) is exact, the returned set of horoballs defines an 
e-ball hull. 

Since (with constant n) the grid G £ is of size Q((s mh.( dx) / s) n -1 ) and computing the convex hull in each flat 



takes 0(M n / 2 J + N log iV) time this proves Theorem 4. 1 



5 Center Points 

In Euclidean space a center point p of a set X C R d of size iV has the property that any halfspace that contains p 
also contains at least N/ (d + 1) points from X. Center points always exist [22 ] and there exists several algorithms 
for computing them exactly [15] and approximately lfT0ll20ll . 

We cannot directly replicate the notion of center points in P(n) with horoballs. Instead we replace it with a slightly 
weaker notion, which is equivalent in Euclidean space. A horo-center point p of a set X C P(n) (or R d ) of size 
N has the property that any horoball that contains more than Nd/{d + 1) points must contain p, where we define 
d = n(n + l)/2 so that P(n) is a d-dimensional manifold. 



Construction for no center point in P(n). Analogous to Euclidean center points, a center point p of X C P(n) 
of N points has the property that any horoball that contains p must also contain at least N/ (d + 1) points from X. 

Theorem 5.1. For a set X C P(n) there may be no center point. 

Proof. Consider a set of distinct points X G P(n) such that all points X lie on a single geodesic a between x\ 
and x jy where x\,x^ G X. Furthermore, let the points lie in a hyperbolic submanifold of P(n). Now, for any 
point p not on a there is a horoball that contains p but contains none of X. So if there is a center point, it must 
lie on q. However, also for any point p G a there is a horoball that intersects a at only p, since the cross-section 
of any horoball in the hyperbolic submanifold will be strictly convex (it can be represented as a hyperball in the 
Poincare model, and geodesies are circular arcs, so there is a horoball tangent to the geodesic at one point). Thus 
for any possible center point p there is a horoball that contains at most 1 point of X. Hence, there can be no center 
point. □ 
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Horo-center points in P(n). The "simple" proof of the existence of center points lfl9l uses Helly's theorem to 
show that a horo-center point always exist, and then in Euclidean space a halfspace separation theorem can be used 
to show that a horo-center point is also a center point. We replicate the first part in P(n), but cannot replicate the 
second part because horoballs do not have the proper separation properties when not defined in R d . 

Theorem 5.2. Any set X C P(n) has a horo-center point. 

Proof. We use the following Helly Theorem on Cartan-Hadamard manifolds (which include P(n)) of dimension 
d ifPTl . For a family of closed convex sets, if any set of d + 1 sets from F contain a common point, then the 
intersection of all sets in 3~ contain a common point. 

In P(n) we consider the family 3~ of closed convex sets defined as follows. A set F G 5F is defined by a horofunc- 
tion b c and a subset X' C X of size greater than Nd/(d + 1) such that X' is the intersection of X and a horoball 
B r (b c ). Then F = B(X') is the ball hull of X', so F C B r (b c ) and F is compact. 

We can argue that any set of d + 1 sets from 3" must intersect. We can count the number of points not in any d + l 
sets as 

d+l d+l 

S < - Nd/(d + 1)) = *£(N(l/(d + 1)) = iV. 

i=l i=l 

So there must be at least one point in X that is in all of the d + l sets. Then by the Helly-type theorem there exists a 
point p such that p G F for any F G 9". 

We can now show that this point p must be a horo-center point. Any horoball that contains more than Nd/ (d + l) 
points from X contains an element of 3~, thus it must also contain p. □ 

5.1 Algorithms for Horo-Center Points 

We provide justification for why it appears difficult to describe an exact algorithm for constructing horocenter points 
in P(n) and then provide an algorithm for an approximate horocenter point in P(n). 

Before we begin we need a useful definition of a family of problems. An LP-type Problem ll24l takes as input a 
set of constraints H and a function uj : 2 H — > R that we seek to minimize, and it has the following two properties. 
Monotonicity: For any F C G C H, u(F) < u(G). Locality: For any F C G C H with uj(F) = uj(G) and 
an h G H such that uj(G U h) > uj(G) implies that to(F U h) > lo(F). A basis for an LP-type problem is a subset 
B C H such that lo(B') < lo(B) for all proper subsets B' of B. And we say that B is a basis for a subset G C H 
if co(B) = to(G) and B is a basis. The cardinality of the largest basis is the combinatorial dimension of the LP-type 
problem. LP-type problems with constant combinatorial dimensions can be solved in time linear in the number of 
constraints 171 . 

Lemma 5.1. A set H of horoballs in P(n), and a function uj(G) = min pg p|^ det(j?) is an LP-type problem with 
constant combinatorial dimension. 

Proof. Monotonicity holds since in adding more horoballs to the a set F C H to get a set G C H (i.e. so F C G) 
we have f| G Q fV 

To show locality, we consider subsets F C G C H such that lu(F) = uj(G). Let P be the set of points 
{p G P(n) | uj(p) = min 9g p| G uj(q)}. Adding a constraint (a horoball) h to G only causes uj(G U h) > u(G) if 

P n h = and thus P n (HguJ = Since He c fV then also p n (fVu/J = $ and u M > 

We now show that our problem has combinatorial dimension d = n(n + 1) /2. P(n) is a d-dimensional manifold. 
Each constraint (a horosphere) is a (d — 1) -dimensional sub-manifold of P(n). Thus letp* = argmin pg p|^ If 
a constraint h lies in the the basis B C F, then p* must lie on the corresponding horosphere. Hence, this reduces the 
problem by 1 dimension. And each subsequent horosphere hi we add to the basis, must also include p*, so it must 
intersect h (and all other horospheres in the basis) trans versally, reducing the dimension by 1. (If h, h' G B do not 
intersect transversally, then we can remove either one from the basis without changing p.) This process can only add 
d horospheres to B because P(n) is d-dimensional, thus the maximum basis size is d. □ 
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This lemma suggests the following algorithm for constructing a horo-center point. Consider all subsets X' C 
X of Nd/(d + 1) points, find the minimal horoball(s) which contain X'. Each of these horoballs can then be 
seen as a constraint for the LP-type problem. Then we solve the LP-type problem, returning a horo-center point. 



Unfortunately, there is no finite bound on the number of horoballs defined by a subset X'. Theorem 3.2 indicates 
that it could be infinite. Thus there are an infinite number of constraints that may need to be considered. 

In order to approximate the horocenter point, we use a similar approach as we did to approximate the ball hull. 
We discretize the set of directions, and create a finite family of constraints for each direction. Then we can solve the 
associated LP-type problem to find a horo-center point. 

More formally, we place a grid G £ on SO(n) so that for any Q' G SO(n) there is another Q G G £ such that the 
angle between Q' and Q is at most e/(^\2y/2sinh(dx/V^))- For each c corresponding to Q G G £ , we consider 
ttf(X), the projection of X onto the (d — l)-flat F corresponding to Q. Within F, we can consider all subsets 
X' C Xof Nd/(d + 1) points and find the hyperplanes defining the convex hull of ttf(X'). This finite set of 
hyperplanes corresponds to a finite set of horoballs which serve as constraints for the LP-type problem in P(n). 

We say a point p is an e-approximate horo-center point if there is a horo-center point p such that for any horo- 
f unction b c we have \b c (p) — b c (p)\ < s. 

Lemma 5.2. A point p that satisfies all of the constraints defined by X and G £ is an e-approximate horo-center 
point of X. 

Proof. Let C(X) C P(n) be the set of horo-center points. Assume that p £ C(X), otherwise let p = p and we are 
done. 

We show the existence of a specific nearby horo-center point p with the following property. Let c be the geodesic 
connecting p and p, and let p = max p / gC <(x) b c (p')- For any q G dC(X) let a be the geodesic connecting q and p. 
If q = argmaxj/gcpn b a (p'), we are done, if not, let q a G dC(X) such that q a = argmax p / gC ( X -) b a (p'). Then 
the geodesic ray on C(X) that goes from q to q a describes a flow on C(X). We can see that the fixed point of this 
flow is p, because as we move to q in the direction of this flow, b a gets closer to is maximum, and the geodesic on 
P(n) from p to q is closer to direction defining the horofunction that q maximizes. 



Now, we can show that b c (p) — b c {p) = 5 < e. This follows by Theorem 4.3 since there must be another direction 
Q G G £ where the corresponding flat contains a geodesic d such that there are more than Ndj (d + 1) points x G X 
such that b(j(x) < b c >(p) and \b c (x) — b c >(x)\ < e. Thus there are more than Nd/(d + 1) points x G X such that 
b c (x) — e < b c (p), and there must be exactly [Nd/ (d + 1) + lj points x G X such that b c (x) < b c (p). Since p and 
p lie on the geodesic c, we have S < e. 



We can now use Lemma 4.6 to bound the difference in values for any horofunction. 1 1 Vb c \ \ is constant for any 
b c , hence for any other horofunction b c / we have \b c i(p) — b c i(p)\ < \b c (p) — b c (p)\ < e (where they are only equal 
when c and d asymptote at opposite points). Since p is a horo-center point, this concludes the proof. □ 



When constructing the set of constraints in each flat F corresponding to a direction in G £ we do not need to 

f N 
\Nd/{dA 



explicitly consider all (atj/zji -n) subsets of X. Each constraint only depends on n points in X, so we can instead 



consider (J) = 0{N n ) subsets of X of size n, and then check if either of the halfspaces it defines in F contain at 
least Nd/ (d + 1) points from X in O(N) time. Only these constraints need to be considered in the definition of p. 

Theorem 5.3. Given a set X C P(n) of size N, (for n constant) we can construct an e-approximate horo-center 
point in time 0((sinh(dx) / 'e) n ~ l N n+l ) time. 

Proof. As per the above construction, for each Q G G £ we only need to consider 0(N n ) potential constraints, 



and each takes O(N) time to evaluate. By Lemma 4.7 G £ is of size 0((sinh(dx)/e) n 1 ) so the total number 



of constraints we need to consider in the LP-type problem is 0((smh(dx) / s) n 1 N n ). The total runtime is thus 
dominated by constructing the constraints and takes 0((sinh(dx) / e)™ -1 N n+1 ) time. □ 
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